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Abstract 

We demonstrate effectiveness of the first-order algorithm from [Mil- 
stein, Tretyakov. Theory Prob. Appl. 47 (2002), 53-68] in application to 
barrier option pricing. The algorithm uses the weak Euler approximation 
far from barriers and a special construction motivated by linear interpola- 
tion of the price near barriers. It is easy to implement and is universal: it 
can be applied to various structures of the contracts including derivatives 
on multi-asset correlated underlyings and can deal with various type of 
barriers. In contrast to the Brownian bridge techniques currently com- 
monly used for pricing barrier options, the algorithm tested here does not 
require knowledge of trigger probabilities nor their estimates. We illus- 
trate this algorithm via pricing a barrier cap let, barrier trigger swap and 
barrier swaption. 

AMS 2000 subject classification. 65C30, 60H30, 91B28, 91B70. 
Keywords. Barrier options, exotic derivatives, weak approximation of 
stochastic differential equations in bounded domains, Monte Carlo tech- 
nique, the Dirichlet problem for parabolic partial differential equations, 
interest rate derivatives. 

1 Introduction 

Barrier option contracts are among the most traded and oldest exotic deriva- 
tives. They accommodate investors' view about the future market behavior 
more closely and they are generally cheaper than the corresponding plain vanilla 
options. Typically, a barrier option is activated (knocked in) or deactivated 
(knocked out) depending on whether a vector of underlying assets or their func- 
tional has crossed a specified barrier level, which itself can be a functional of the 
underlying assets. Due to its attractive features, barrier optionality has been 
introduced in a wide range of derivatives products. In the context of credit 
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risk, the event of default of some reference entity can be modelled as a lower 
barrier on the equity of the entity. This is the key idea in structural models 
for pricing popular credit instruments such as credit default swaps (CDSs) and 
credit default obligations (CDOs). More recently, barrier optionality has been 
also used to model contingent convertible (CoCo) bonds which were introduced 
to provide financial institutions with sufficient capital in times of distress or 
systemic risk. 

Closed form solutions for barrier option prices can be obtained only in some 
particular settings. For instance, they are available in the case of a single un- 
derlying asset and a constant barrier within the standard Black-Scholes setup 
(see, e.g. [TU] HH1 HI E3 121] )• in products involving a large number of depen- 
dent assets numerical approximation for pricing and hedging barrier options is 
usually inevitable and this can be a challenging problem. 

In this paper we assume that underlying assets are modelled via multidimen- 
sional stochastic differential equations (SDEs) and we consider European-type 
barrier options. The arbitrage price u(t, x) of such an option solves the Dirichlet 
problem for a linear parabolic partial differential equation. Finding this price 
numerically requires efficient weak approximations of diffusions in a bounded 
domain. 

"Ordinary" numerical methods for SDEs on a finite time interval in R d 
are based on a time discretization HH [IB] . They ensure smallness of time 
increments at each step, but might not ensure smallness of space increments. 
In [T3l [T2l [15] (see also [16]) a number of weak-sense approximations for SDEs 
in a bounded domain were proposed, in which space increments are controlled 
at each step so that the constructed approximation belongs to the bounded 
domain. Approximations of pj] (see also [16]) are based on adaptive control of 
a time step of numerical integration of the SDEs. A step is chosen such that (of 
course, aside of reaching a required accuracy) the next state of a Markov chain 
approximating in the weak sense the SDEs' solution remains in the bounded 
domain with probability one. This leads to a decrease of the time step when 
the chain is close to the boundary of the domain. The chain is stopped in a 
narrow zone near the boundary so that values of the solution u(t,x) (i.e., the 
option price at time t and underlyings' price x) in this zone can be approximated 
accurately by the known values of the function tp on the boundary (i.e., the value 
of the option at the barrier). Another type of approximations was proposed in 
[15] (see also [16]). In the algorithm of weak order one from [15] the step of 
numerical integration of the SDEs is constant for points belonging to a certain 
time layer t = tk- Far from the boundary, a Markov chain approximating the 
SDEs' solution is constructed using the weak Euler scheme (i.e., using discrete 
random variables for approximating the Wiener increments). When a point is 
close to the boundary, we make an intermediate (auxiliary) step of the random 
walk, which preserves the point in the time layer t — t^. On this auxiliary step 
we "flip a coin" to decide whether to terminate the chain on the boundary or 
jump back in the domain and continue the random walk. The construction of 
this step is based on the idea of linear interpolation for the solution u(t, x). The 
algorithm is efficient and very easy to implement. Its simpler version of order 
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1/2 in the weak sense is also presented in [I3J[TB]. In Section [5] we recall these 
two algorithms from [T5l [16] . Due to our knowledge, despite simplicity of these 
weak schemes, they have not been used in financial applications. In this paper 
we try to fill this gap and illustrate their applicability to pricing barrier options. 

Currently, the popular numerical approach for pricing barrier options ex- 
ploits the Brownian bridge technique [TJ [H [501 (also see [S] for a review and 
the references therein) . It is based on simulation of a one-dimensional Brownian 
bridge extremum between time steps and computing analytically the associated 
probability of exiting the spatial domain for each time interval of the partition. 
It was proved in [7] that this approach realized along with the Euler scheme 
(which uses Gaussian random variables for simulating Wiener increments) re- 
sults in an approximation of weak order one. The Brownian bridge technique 
relies on analytical formulas for trigger probabilities and can run into difficulties 
in the case of multiple barriers and/or correlated structure of the underlyings 
when there are no closed formulas for the distribution of extremum. Though 
some extensions to these more general and not uncommon problems have been 
considered, e.g. in [51 |2"U]. 

In contrast the simplest random walk algorithm of [15] displays a high degree 
of flexibility and can be applied to various structures of the contracts including 
derivatives on multi-asset correlated underlyings and can deal with various type 
of barriers, e.g. single, double and time dependent barriers. In comparison 
with the Brownian bridge techniques the method of |15j does not require the 
knowledge of the trigger probabilities nor their estimates. 

In Sections 0][6] we present three examples on how to apply the algorithm 
from [15] for valuation of barrier options. These contracts cover the most com- 
mon types of barrier options. In the first example (Section [4]), we consider 
an algorithm for barrier derivatives where the payoff depends on a single un- 
derlying. As an illustration, we deal with pricing a barrier caplet and provide 
ready-for-implementation procedure which can be easily applied to similar other 
contracts. The second example (Section[S]) is devoted to multi-asset options with 
barriers imposed on all or some of the correlated underlying assets. We illus- 
trate this case by pricing a trigger swap. The last example (Section [5]) is barrier 
contracts written on an asset that can be expressed through some other multi- 
asset underlying. As a specific case, we consider valuation of a barrier swaption 
under the LIBOR market model (LMM). 

2 Simplest random walks for stopped diffusions 

Let F, P) be a complete probability space, Ft, < t < T, be a filtration 
satisfying the usual hypotheses, (wt, Ft) be an r-dimensional standard Wiener 
process. Let G be a bounded domain in R d and Q = [to,T)xGbe a cylinder in 
R d+1 , L = Q\Q be the part of the cylinder's boundary consisting of the upper 
base and lateral surface. Price of barrier options with underlying modelled by 
a diffusion process can usually be expressed as 

u(t,x) = E[(p(T,X t , x (r))Yt, x ,i(T) + Z t , x ,i,o(r)] , (2.1) 
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where X tiX (s), Y t . x . y (s), Z tiXi y iZ (s), s > t, is the solution of the Cauchy problem 
for the system of SDEs: 

dX = (b(s,X)-a(s,X)fi{s,X))ds + a{s,X)dw(s), X(t)=x, (2.2) 
dY = c(s,X)Yds + fj, T (s,X)Ydw(s), Y(t)=y, (2.3) 
dZ = g(s,X)Yds + F T (s,X)Ydw(s), Z(t) = z, (2.4) 

(t, x) £ Q, and r = Tt, x is the first exit time of the trajectory (s, X ttX (s)) to the 
boundary T. In (|2.2p - (12.4[) . b(s, x) is a d-dimensional column- vector, the a(s,x) 
is a d x r matrix, /i(s, x) and -F(s, a; ) are ^-dimensional vectors, and Y(s), Z(s), 
c(s, X) and g(s, X) are scalars. We assume that all the coefficients in (|2.2I) - (|2.4[) , 
the function ip(t, x) defined on T and the boundary dG of the space domain G 
satisfy some regularity conditions. 

We note that the value of the expectation u(t, x) in (|2.1[) does not depend 
on a choice of functions fi(s,x) and F(s,x). This flexibility can be used for 
reducing variance of the random variable ip(r, X t , x {T))Y t , x ,i{T) + Z tjX ^.o(T) with 
the aim of reducing the statistical error in computing u{t, x) via the Monte Carlo 
technique [TTJ [TB] . For instance, if [i and F are such that 

E ^ + V + F j =0, j = 1, . . . , r, (2.5) 

then Var[(p(T,Xt )X (T))Y tlXt i(T) + Zt.^i.oW] = and ^(r,Xt, x (T))Ft, x ,i(r) + 
■^t,a;,i,o(' r ) = u(t,x) [HI [T^. As we see from (I2.5[) . optimal \i and F require 
knowledge of the the solution u(t,x) (i.e., the option price) to the considered 
problem and its derivatives (i.e., deltas) which is impractical. However, instead 
of the exact u(t, x) in (|2.5p . one can use its approximation (e.g. price and deltas 
for a related option for which the closed-form solution is known) to find some 
suboptimal /i and F which can lead to variance reduction [SJ 117) . 

To simulate (|2.ip - (|2.4[) , we need an approximation of the trajectory (s, X{s)) 
which satisfies some restrictions related to its nonexit from the domain Q. Let 
us recall two algorithms for (|2~T|) - ([2T4l) from [T51H6]. 

We apply the weak explicit Euler approximation with the simplest simulation 
of noise to the system (f2"T2"T) -(f2~4 |) : 

X t . x {t + h) « X = x + h(b(t,x)-a(t,x)n(t,x)) + h 1/2 a(t,x)t, (2.6) 
Y t , x , y (t + h) w Y = y + hc(t,x)y + h^ 2 ^(t,x)y^ , (2.7) 
Z t:X:V:Z {t + h) a Z = z + hg(t,x)y + h 1 / 2 FT(t,x)yt, (2.8) 

where h > is a time-discretization step (a sufficiently small number), £ = 
(£ , . . . , £ r ) T , C i i — 1; • • • j r i are mutually independent random variables taking 
the values ±1 with probability 1/2. Clearly, the random vector X takes 2 r 
different values. 

Introduce the set of points close to the boundary (a boundary zone) St,h C G 
on the layer t : we say that x € S 1 ^^ if at least one of the 2 r values of the vector 
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X is outside G. It is not difficult to see that due to compactness of Q there is a 
constant A > such that if the distance from x £ G to the boundary dG is equal 
to or greater than X\/h then x is outside the boundary zone and, therefore, for 
such x all the realizations of the random variable X belong to G. 

Since we should impose restrictions on an approximation of the system (|2.2I) 
so that it does exit from the domain G, the formulas (|2.6[) - (I2.8[) can be used 
only for the points x £ G\S t ,h on the layer t, and a special construction is 
required for points from the boundary zone. Let x £ S t .h- Denote by x* £ dG 
the projection of the point x on the boundary of the domain G (the projection 
is unique assuming that h is sufficiently small and dG is smooth) and by n(a: 7r ) 
the unit vector of internal normal to dG at x n . Introduce the random vector 
X£ h taking two values x v and x + h 1 ^ 2 Xn(x 7r ) with probabilities p — p x ,h and 
q = q x ,h = 1 - Px,h, respectively, where 

Ph - ^ (29) 

Px ' h ~ \x + hV 2 \n(x*)-x*\ ' K ' 

This construction is motivated by the following observation |15j . If v{x) is a 
twice continuously differentiable function with the domain of definition G, then 
an approximation of v(x) by the expectation Ev(X™ h ) corresponds to linear 
interpolation and 

v(x) = Ev(X* th ) + 0(h) = pv(x*) + qv(x + h 1/2 \n(x n )) + 0(h) . (2.10) 

We emphasize that the second value x + h 1 ^ 2 \n(x^) does not belong to 
the boundary zone. We also note that p is always greater than 1/2 (since the 
distance from x to dG is less than /i 1 / 2 A) and that if x £ dG then p = 1 (since 
in this case x n = x). 

Let a point (to, Xg) £ Q. We would like to find the value u(to, Xq). Introduce 
a discretization of the interval [^ch^l > f° r definiteness the equidistant one: 

t < h < ■ ■ ■ < t M = T, h:=(T- t )/M. 

To approximate the solution of the system (|2.2p . we construct a Markov 
chain (tk,Xk) which stops when it reaches the boundary L at a random step 
yc < M. The resulting algorithm can be formulated as Algorithm 12.11 given 
below. 

It is proved in [15] (see also [16]) that under appropriate regularity assump- 
tions on the coefficients of (|2.2p - ()2.4|) . the boundary condition ip(t,x) in (|2.1[) 
and on the boundary dG Algorithm 12.11 converges with weak order one. 

The next algorithm is obtained by a simplification of Algorithm 12. II as soon 
as Xk gets into the boundary domain St k ,h, the random walk terminates, i.e., 
7t = k, and X^ = X£, Y M = Yj,, Z M = is taken as the final state of the 
Markov chain. The resulting algorithm takes the form of Algorithm 12.21 

It is proved in j!51 116] that under appropriate regularity assumptions on the 
coefficients of (|2.2j) - ([2.4[) , the boundary condition (p(t,x) in (|2.1[) and on the 
boundary dG Algorithm 12.21 converges with weak order 1/2. We note that in 
one-dimension (i.e., in the case of a single underlying) Algorithm l2.2l is analogous 
to pricing barrier options by binary trees (see, e.g. [5]). 
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Algorithm 2.1 Algorithm of weak order one for (|2.1[) - (I2.4I) 
STEP 0. X' = x Q , Y = 1, Z = 0, k = 0. 

STEP 1. If X' k i S tk ,h, then X k = X' k and go to STEP 3. 

If X' k £ St k ,h, then either Xk = X k with probability 
p X ' k ,h or X k = X' k + h lI2 \n{X'^) with probability q x > k ,h ■ 

STEP 2. If X k = Xg, then STOP and x = k, 

X^ = X£ , = Yk , Z^ = Zk ■ 

STEP 3. Simulate^ and find X k+1 , Yfc+i, Zk+i according to 

for t = t k , x = X k , y = Y k , z = Z k , 

STEP 4. If k + 1 = M, STOP and x = M, X„ = X' M , Y„ = Y M , 
Z^ = Zm, otherwise k = k + 1 and return to STEP 1. 



Algorithm 


2.2 Algorithm of weak order 1/2 for (j2Tj)-(j2~4)) 


STEP 0. 


X = xo, Y = 1, Zo = 0, k = 0. 


STEP 1. 


If X k i S tk , h , then go to STEP 2. 




If X k £ S tk , h , then STOP and x = k, X„ = XI, 




Yft = Yk, Z& = Zk ■ 


STEP 2. 


Simulate £ fc and find X k +i, ifc+i, Zk+i according to 




$ZM-$n$ for * = t k , x = X k , y = Y k , 




z = z k, € = £fc • 


STEP 3. 


If k + 1 = M, STOP and x = M, X^ = A M , ^ = Ym, 




= Zm, otherwise k — k + 1 and return to STEP 1. 



3 LIBOR Market Model 

We will now assume that there exists an arbitrage-free market with continuous 
and frictionless trading taking place inside a finite time horizon [i ,t*]. 

Among the most important benchmark interest rates is the London Inter- 
bank Offered Rate (LIBOR). It is based on simple (or simply compounded) 
interest. The forward LIBOR rate L(t, T,T + S) is the rate set at time t for the 
interval [T,T + 5], t < T. If we enter into a contract at time t to borrow one 
unit at time T and repay it with interest at time T + 5, the interest due will be 
5L(t,T,T + S). 

A simple replication argument (see, e.g., |4]) relates LIBOR rates and bond 
prices via the following identity 

where Pit, T) is the price at time t < T of a default-free zero coupon bond. 
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For simplicity, we fix an equidistant finite set of maturities or tenor dates 



T <- - <T N = T*, Ti = i6, i = 0,...,N, 



(3.2) 



where 



5 = (T* — T )/N, 



denotes the fixed length of the interval between tenor dates. 

Let us introduce a simplified notation for the time t forward LIBOR rate for 
the accrual period [7i,Ti_|_i] and the payment at Tj+i: 

L\t) : = L(t,Ti,T i+1 ), 

to < i<fAT„ h<T t <T*, i = 0,...,N-l. 

In the case of LIBOR Market Model (LMM) the arbitrage-free dynamics of 
L l (t) under the forward measure Q T * +1 associated with the numeraire P(t, T^+i) 
can be written as the following system of SDEs (see, e.g. [H [TBI |2"T]): 



dL l {t) 



E iSzSfepij^W* + ^(t)dw? k+1 (t), t >k 7 t<T k , 
j=k+l 



= < 



<n{t)dWi k+1 {t), i = k, t<T h 



(3 ' 3) 

where W Tk+1 = (W Q Tk+1 , . . . , W^ k ^) T is an iV-dimensional correlated Wiener 
process defined on a filtered probability space ^f2, J 7 , {J~t} ta<t<t , , Q Tfc+1 ); the 
instantaneous correlation structure is defined as 



rT k 4 



E 



II 



-(*)Wf fc+1 (t) 



Pi,,-, *>i = o, 



,iV-l: 



(3.4) 



and cri(t), i = 0, . . . , N — 1, are instantaneous volatilities which we assume here 
to be deterministic bounded functions. 

Let p be the instantaneous correlation matrix with elements p i j . To simulate 
the correlated Wiener processes, we will use the pseudo-root of the correlation 
matrix p defined via the equation 



p=UU\ 



(3.5) 



where U is an upper triangular matrix. Using U and introducing iV-dimensional 
standard Wiener process, one can re- write (|3.3|) in the form of (|2.2|) . 

In what follows we will assume that the current time to is set to 0. For con- 
venience, we also assume a unit notional value of all the contracts we introduce 
below. In our numerical experiments in the next sections we take the correlation 
function of the form: 

Pid =ex9(-P\T i -T J \). (3.6) 
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4 Barrier cap/floor 



In this section we consider Monte Carlo evaluation of barrier options written 
on a single underlying. We use a knock-out caplet for illustration, though our 
treatment is rather general and can be used to value different barrier option, 
for instance European and Parisian barrier options and options with different 
barriers including single, double and time-dependent barriers, both for fixed- 
income and equity markets. 

An Interest Rate Cap is a security that allows its holder to benefit from 
low floating rates and be protected from high ones. Similarly, an Interest Rate 
Floor is an instrument designed to protect from low floating interest rates yet 
allow the holder to benefit from high rates. Formally, a cap price is obtained 
by summing up the prices of the underlying caplets, call options on successive 
LIBOR rates. Also, a floor is a strip of floorlets, put options on successive 
LIBOR rates. 

A knock-out caplet pays the same payoff as a regular caplet as long as a 
prescribed barrier rate H is not reached from below by the corresponding LIBOR 
rate before the option expires. More specifically, the price at time t < Tq of 
knock-out caplet set at time Ti_i with payment date at Ti, i > 1, with strike 
K and unit cap nominal value is given by 

V cap i et (t) = SP(t,T l+1 )E QT ' +1 [(L*(r<) - K) +x (8 > T<)| Ji] , (4.1) 

where 6 is the first exit time of L' l (s), s > t, from the interval G = (0,H). Let 
t be the first exit time of the space-time diffusion (s,L l (s)) from the domain 
Q = [t, Ti) x (0, H). Obviously, t = 9AT,. 

The dynamics of L l (s) under Q Ti+1 is (see (|3.3jl ): 



dL l 



ai{s)dW? i+1 (s), s<T t . (4.2) 



Note that the correlation structure of (|3 . 3|) does not influence the price of the 
knock-out caplet since it does not depend on the joint dynamics of forward rates. 

One can observe that the dynamics (|4.2p coincides with the model of a stock 
price process under the risk-neutral measure in the case of zero interest rate. 
This means that by dropping the factor 6P(t, Tj+i) in (|4.ip . the valuation of 
European barrier options on equity with zero interest rate and of barrier caplets 
under the LMM coincide. 

In the considered case the price of the barrier caplet has the well-known 
closed-form solution: 

V ca plet(t) = V cap l et {t,L l {t)) (4.3) 

= 5P(t,T l+1 ) {L*(t) [<S>(5 + (V(t)/K,v t )) - <f>(5 + (V(t)/H, Vl ))] 
-K K, Vi)) - *(*_(£*(*)/#>«<))] 

-H [^(S + (H 2 /(KL\t)),v t )) - «0)] 
+KL\t) [<S>{5-(H 2 /{KUit))^)) - <5>{5_{H/L\t), Vi ))\ /H) , 
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where $(•) denotes the standard normal cumulative distribution function and 
5+(x,v) = (lnx + v 2 /2)/ Vl 5_{x,v) = Qnx-v 2 /2)/v, (4.4) 

and 

vf = J (<7i(s)) 2 ds. 

This analytical result will be used in our numerical experiments to access the 
performance of proposed algorithms. We note that the algorithm presented 
in this example can be easily extended to a more general model of underlying 
when the closed-form solution might be not available. In particular, there is no 
difficulty in including a drift term in the underlying dynamics (see also Sections[5] 
and [5]). In the experiments we simulate 

V caplet (t) = V cap i et (t)/SP{t,T l+1 ), (4.5) 

i.e. we drop SP(t,Ti + i) from (|4.1[) . which does not imply any loss of generality 
since the caplet price can easily be recovered by multiplying V cap i e t(t) by the 
factor 8P(t,Ti + i) observable at time t. 

4.1 Algorithm 

To preserve positivity of the LIBOR rate, we simulate the log dynamics cor- 
responding to (|4.2p rather than the LIBOR rate L l {t) itself. To illustrate the 
variance reduction technique discussed in Section [5J we complement (14.21) with 
the equation (cf. I|2.4p): 

dZ = F{s,L l )dW^ +1 {s), Z(0) = 0, (4.6) 

with (see fl23J and 

F( Sl V) = -a t ( S )^-V caplet (s,L l ( S )) . (4.7) 

We choose a time step ft. > so that M = Tijh is an integer. We set InLg = 
L l (0) and Z Q = 0. The weak Euler scheme (EU, {2^) applied to ((H) in the 
log form and (|4.6p takes the form: 

\nL l k+1 = InLl-^a^fh + a^Vh^^, (4.8) 

Zk+i = Z k + F( S ,Li)Vh£ k+1 , (4.9) 

where £ fc are independent random variables distributed by the law P(£ = ±1) = 
1/2. 

The boundary zone St : h required for Algorithms 12.11 and 12.21 is chosen here 

as 

St hth = {L\: InLl^lnH+^aUt^h-a^Vh}, (4.10) 
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i.e., the condition for hiL k+1 to be inside the domain G is 



In4 < In if + -of (tjfc)h - (n(tk)Vh , (4.11) 

and the corresponding Afe for Algorithm 12.11 is so that 

\ k Vh = -^(t k )h + a l {t k )Vh. . (4.12) 

We note that instead of (|4.10[) and (|4.12[) we could take the wider boundary 
zone St k ji = {L l k : lnL' fc > \nH — <Ti(t k )Vh} and correspondingly \ k = Ui(t k ). 
A wider boundary zone usually leads to a bigger numerical integration error. 
In this example we cannot take a boundary zone narrower than St k ,h in (|4. 10[) 
because it would not ensure that the chain In L k belongs to G. 

To realize Algorithm 12. 1[ we follow the random walk generated by (|4.8j) and 
at each time t k , we check whether at the next step L k+1 cannot cross the barrier 
if, i.e., we check whether the condition (j4.11j) holds. If it does, we perform 
(|4.8[) - (|4.9I) to find ]n.L\ +1 , Z k+i . Otherwise, L\ has reached the boundary zone 
St kt h, where we make the auxiliary step: we either stop the chain at In if with 
probability p : 

X k Vh 

P ~ lnF-lni^ + AfcV^ 

or we kick the current position of the random walk In L k back into the domain 
to the position In L\ — \ k V~h with probability 1 — p and then carry out (|4.8|) - (|4.9p 
to find hiL k+1 , Z k+ \. If k + 1 = Ai, we stop, otherwise we continue with the 
algorithm. The outcome of simulating each trajectory is a point (t^, In L^, Z^). 

To realize Algorithm l2.2[ we also follow the random walk generated by (|4.8j) . 
and at each time t k , we check whether the condition f|4. 1 1[) holds. If it does not, 
L l k has reached the boundary zone St kl h and we stop the chain at In H. If it does, 
we perform (|4.8[) - (|4.9[) to find \nL l k+1 , Z k+1 . If k + 1 = M, we stop, otherwise 
we continue with the algorithm. The outcome of simulating each trajectory is 
again a point (t^, hi-L^., Z^). 

In the experiments we evaluate the expectation 

V cap i et (0) = £ qT * +1 [(L\T 1 )-K) +x (6>T 1 ) + Z(t)] (4.13) 

w £ QT * +1 [(exp(ln Li) -K) +X (x= M) + Z» 

The approximate equality in (|4.13|) is related to the bias due to the numerical 
approximation. The expectation on the right-hand side is realized via the Monte 
Carlo technique. 



4.2 Numerical results 

Here we present some results of numerical tests of Algorithms 12.11 and 12.21 for 
pricing the barrier caplet (|4.13[) . We use the following parameters in the exper- 
iments: i = 9, K = 1%, H = 28%, (5 = 1, L 9 (0) = 13%. The volatility a l {t) is 
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assumed to be constant at 25% . The exact caplet price with these parameters 
evaluated by (|4.3[) is 6.57%. In the experiments we did 10 6 Monte Carlo runs. 
The results are presented in Figure 14.11 We see that Algorithm 12.11 is much 
more accurate than Algorithm 12.21 We also observe "oscillating" convergence 
which is typical for binary tree methods [S]. 



Figure 4.1: Barrier caplet price: Comparision of the results of numerical exper- 
iments for the Algorithm 12.11 (Algorithm 0(h)) and Algorithm 12.21 (Algorithm 
0(Vh)) and the exact caplet price (solid line) evaluated for i = 9, K = 1%, 
H = 28%, 5 = 1, L 9 {0) = 13%, <7 % {t) = 25%. 

Let us also remark on the effect of variance reduction in these experiments. 
For instance, in Algorithm 12.11 for h = 0.02 we got the Monte Carlo error 
(i.e., half of the size of the confidence interval for corresponding estimator with 
probability 0.95) equal to 1.11 x 10~ 4 in the case of F = and 1.55 x 10~ 5 in 
the case of the optimal F from (|4.7[) (i.e., 100 time speed-up in reaching the 
same level of the Monte Carlo error). The use of the optimal F does not result 
in zero Monte Carlo error due to the error of numerical integration. 
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5 Trigger swap 



This example is devoted to evaluation of multi-asset barrier options with barriers 
on all or some of the correlated underlying assets. We consider a trigger swap 
as a specific case, though the considered approach can be used to value other 
multi-asset barrier options, for instance basket options, CDOs and n t,l -to-dcfault 
CDSs, and it can also be applied for options with single, double and time- 
dependent barriers. 

A trigger swap is a swap on a floating reference rate that takes effect or 
terminates when some index rate hits a specified trigger level. Trigger swaps 
have a number of variations [H [3T] . Here we consider a knock- in version of a 
payer trigger swap with a fixed rate K whose barrier is continuously monitored. 
The index and reference rate both coincide with a LIBOR rate. For given trigger 
levels H ,..., H 1 *- 1 associated with the LIBOR rates L°(t), L N - X (t), the 
structure of the swap under consideration is expressed as follows. Once one 
of the the continuously monitored LIBOR rate L°(t), . . . , L A ' _1 (t) for the first 
time hits the corresponding trigger level H°, . . . , H 1 ^^ 1 from below, the contract 
holder enters into the payer swap starting at next tenor date for the remaining 
time to the last tenor Tn-i- More specifically, let 8 be the first exit time of 
L°(s), L^-^s), s > 0, from the domain G = (0, H°) x • • • x (0, iJ^" 1 ) , 
r be the first exit time of the space-time diffusion (s, L°(s), . . . , L N ^ 1 (s)) from 
the domain Q = [0, Tn-i) x G (clearly r = 9 A Tjv-i), and T e ( T ) be the closest 
tenor date Ti to r from the right, i.e., g(t) is defined as 



g(t) 



i{i, i = 0,l,...,N-l:t<Ti}. 



If 6 < Tjv_i, then at a tenor date the contract holder enters into the 

contract according to which the holder pays fixed payments of 5K and receives 
floating payments of <5L ,-1 (2i_i) at the coupon dates Tj, i = q(t) + 1, . . . ,N; 
otherwise the contract expires worthless. 

The value of this trigger swap at time t = under the forward measure Q TjY 
is given by 



Vi 



trswap 



(0) 



p(o, t n )e qT * 



1 



P{T g ( T ),T N ) 
l-P(T e{T) ,T N )-KS 

or in terms of the LIBOR rates 



(5.1) 



N 

E 

»=8(t)+1 



P(T e{T) ,Ti) X (0<Tn-i) 



Vtrswap 

(0) = F(0,T w )^Q A 



N-l 



N N-l 



(5.2) 



- K6 E II i 1 + SLj (T e{ r))) - 1 J X(6 < Tn-!) 

>=s(t)+1 j=i 
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In order to price this contract by the Monte Carlo technique, we need to 
generate paths for the vector L e ^(t), . . . , L N ~ x (t). This means that the size of 
the vector of the LIBOR rates which we need to simulate decreases over time. 
The dynamics of L l (t) under Q Tjv are described by the SDEs (cf. 



^- = - ai (t) i S +8lJ(t) *W (*)<ft+gi(* (t), i = g(t),...,N-l. 

(5.3) 



5.1 Algorithm 

Here we only apply Algorithm 12.11 For simplicity, we consider such set of tenor 
dates Ti and time steps h that Ti/h arc integers. 

As before, we simulate dynamics of the LIBOR rates L e ^(t), . . . ,i JV_1 (t) 
in log space according to the weak Euler scheme (cf. ()2.6[) ): 

1 + SL J k 

1 2 r N 

Z 3 = Q k 

i = Q k+ i,--',N -1, 

where £j k are mutually independent random variables distributed by the law 
P{£ = ±1) = 1/2 and g k =: g(t k ). 

The algorithm for the considered trigger swap proceeds as follows. Let L k = 
(L e k k , . . . , L^" 1 ) 1 . Denote by >c the first exit time of (t k , L k ) from Q. Let M — 
Tff-i/h. Suppose by time step k none of the rates L e k k , . . . ,L k _1 have crossed 
their barriers H ak , . . . , , i.e., \ { x ^ ^) — 0- Then we evaluate whether 

at the next time step k + 1 the event yc = k + 1 might be realized. One can see 
that the rate L l k+1 , i = g k+1 , . ■ ■ ,N — 1, computed via (|5.4p will be below the 
barrier H % , i.e. inside the domain G, if the following is true 

\nL\ < lniT - X k Vh , (5.5) 

where 

Afc = CT Max \/N — Q k 

and a Max = maxj^ <Jj{t k ). 

If ()5.5[) is satisfied at time t k for all rates L k k+1 , . . . , L k ~ x then we move to 
step tfc+i, evaluate kiL k+1 , i — g k+1 , . . ■ , N — 1, according to (|5.4[) and continue 
with the algorithm unless k + 1 = M (in this case the trigger swap expires 
worthless). 

The case when the condition (|5.5j) does not hold for a single i implies that 
the point L k is in the boundary zone St kt h and is near the barrier H l . Then we 
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either assign x = k, lni^ = \nH l , hiL^ — \nL k for j ^ i, T e ^ = T 0k with 
probability 

p = : : t= (5.6) 

In if 1 -lni| + A fcX /7T 



and carry on with simulating lnLJ, +1 , i = Q^+i, . . . ,N — 1, according to (|5.4[) 
starting from In until time T g ^ = min {T^ : < Tj, i = 0, 1, . . . , N — 1} (the 
barriers are "removed" in simulating the remaining part of this trajectory); or 
we jump outside the boundary zone St kl h by changing the i th component of 
In Lfc from In L l k to In L l k — XkVh with probability 1 — p l , perform the usual step 
according to (15.4[) and continue with the algorithm unless k + 1 = M (in this 
case the trigger swap expires worthless). 

We note that in comparison with the original formulation of Algorithm 12.11 
here we do not stop the chain Lk at its first exit time from the space domain 
G. Instead, when the barrier is hit, we find the trigger tenor date T e , and if 

< Tjv_i, we continue the simulation according to (|5.4I) until T Q to get the 
required LIBOR rates V(T e{T) ) in 

Now, we discuss the case when the condition (|5.5|) does not hold for more 
than one i (i.e., the random walk has reached a corner of the domain G). In this 
case, the algorithm proceeds as follows. Let us denote by I = {l±, . . . , l n } the 
set of tenor dates corresponding to the LIBOR rates for which (|5.5|) is violated. 

First, we select the rate from the set jlnl^. 1 , . . . , InLj™ j , which is the closest 

to its boundary, i.e. lj such that \nH lj — lnL^ is minimum over j = 1, . . . ,n. 
Then, we repeat the procedure which is given above for the single i with the 
following difference. If hiL^ jumps from the boundary to haLl — X^Vh, we find 

the second closest rate from the set |ln L l k , . . . , ln l) k j and as before repeat 
for this point the routine we have presented for the single i. We follow this 
procedure in the outlined fashion until either the set £ is empty or for some lj 
we reach the boundary and assign 

The outcome of simulating each trajectory is the payer swap starting tenor 
date T e , the stopping time h and the point lnL^ with r\ — T e>c /h, which are 
used for evaluating the trigger swap: 



V tr swap 

(0) « P(0,T N )EV 



N 



N-l 



n (i + si*) 



N N-l \ 

~ K6 E II ( x + SL i) - 1 *(* < M ) 

with the expectation simulated by the Monte Carlo technique. We present 
the pseudocode for simulating a single trajectory based on the algorithm we 
described above. 
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Algorithm 5.1 Pseudocode for simulating a single trajectory in pricing the 
barrier trigger swap 
SET M to Tjv_i//i, k to 1, k to M 
WHILE k < M 
IF k > k 

FOR j = to N - 1, 

IF (|5.5[) for j is false 

calculate p-* by (|5.6p 
form array p of p- 7 
ENDIF 
ENDFOR 

sort p in descending order 
FOR n = 1 to length(p) 

generate u ~ Unif [0.1]; 
IF u < p(n) 

SET k to A: 
SET InL| to lniT 
SET M to T Bk /h 
BREAK 
ELSE 

SET In 14 to hiLi-XkVh 
ENDIF 
ENDFOR 
ENDIF 

Evaluate lnL fe+ i by (|5.4p 
Increase k by 1 
ENDWHILE 



5.2 Numerical results 

Let us present results of numerical experiments we performed for pricing a 
trigger swap with Algorithm 16. 11 The parameters chosen for the experiments 
are T° = 5, T* = 16, K = 0.01, H = 0.13, 5 = 1, /3 = 0.2. The initial LIBOR 
rate curve is assumed to be fiat at 4% and the volatility (Tj(t) is set to be constant 
at 20% . In the simulations we run 10 6 Monte Carlo paths. 

Since the closed-form formula for trigger swap (|5.2[) is not available, we found 
the reference trigger swap price by evaluating the price using Algorithm 15.11 
with h = 0.01 and the number of Monte Carlo runs 10 6 . This reference price is 
5.46 x 10" 2 with the Monte Carlo error 5.50 x 10" 4 , which gives half of the size 
of the confidence interval for the corresponding estimator with probability 0.95. 

The results of the experiments with Algorithm l5 . H are presented in Table I5TT1 
In the table, the values before "±" are estimates of the bias computed as the 
difference between the reference price and its sampled approximation, while 
the values after "±" give half of the size of the confidence interval for the 
corresponding estimator with probability 0.95. The "mean exit time" is the 
average time for trajectories (tk,Lk) to leave the space-time domain Q. The 
experimentally observed convergence rate for Algorithm lS.ll is in agreement with 
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the theoretical first order convergence in h (though we note that the convergence 
theorem in [15l [16] is proved under restrictive regularity conditions and the 
payoff of the trigger swap and the boundary of the space domain G do not 
satisfy these conditions). 



Table 5.1: Performance of Algorithm 5.1 for the trigger swap. 



h 




error 




mean exit time 


0.25 


2.22 x 10 


~ 2 ± 6.39 x 10" 


-4 


12.51 


0.2 


1.85 x 10 


-* ± 6.26 x 10" 


-4 


12.61 


0.125 


1.17 x 10 


-* ± 6.01 x 10" 


-4 


12.78 


0.1 


9.56 x 10 


~' A ± 5.92 x 10" 


-4 


12.83 


0.0625 


6.03 x 10 


~' A ± 5.78 x 10" 


-4 


12.92 


0.05 


4.67 x 10 


- 3 ±5.72 x 10" 


-4 


12.95 



6 Barrier swaption 

In this section we consider Monte Carlo evaluation of a knock-out swaption 
under the LMM. We use the knock-out swaption as a guide in our exposition, 
its treatment is rather general and it can be used to value different barrier 
options, where the underlying and barrier can be expressed as functionals of 
some diffusion process. 

A European payer (receiver) swaption is an option that gives its holder a 
right, but not an obligation, to enter a payer (receiver) swap at a future date 
at a given fixed rate K . Usually, the swaption maturity coincides with the first 
reset date To of the underlying swap. The underlying swap length Tv — To is 
called the tenor of the swaption. 

Without loss of generality, we concentrate on a knock-out receiver swaption 
with the first reset date Tq. A knock-out swaption has the structure as a stan- 
dard swaption except that if the underlying swap rate is above a barrier level 
R up at any time before Tq then the swaption expires worthless. The price of the 
knock-out swaption at time t — under the forward measure Q T ° is given by: 



^swaption 

(0) = p(o,t )eQ To 



N 

S (R swap (T ) - K) + J2 P(T , Tj) X (0 > T ) 
j=i 

(6-1) 

where 9 is the first exit time of the process R SW ap{s), s > 0, from the interval 
(0,R up ). The swap rate R SW ap{s) can be expressed in terms of the spanning 
LIBOR rates as 

N-l 

l-v n (i+^'oo) 

Rs Wap (s) = JV _ 1 . (6.2) 

*£vn (i+*#oo) 

l= o i=° 
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The bond prices P(T ,Tj) can also be expressed via LIBOR rates (see (|3.ip ). 

Also, let r be the first exit time of the space-time process (s, R S wap(s)) from 
the domain D = [0, T ) x (0, R up ) (obviously, r = 6 A T ). 

We note that expression (|6.1|) depends on the joint distribution of the forward 
rates L°(T ), L^-^To). The LMM dynamics of LIBOR rates under Q T ° are 
given by (cf. (1331) ): 

= "•<(*) E iSft) ^^^^ + i = 0,...,N-l. 

(6.3) 

In this example we deal with pricing the barrier swaption (|6.ip expressed in 
terms of the spanning LIBOR rates with dynamics in the form of (|6.3I) . This 
means that we consider this problem in the coordinate system of the spanning 
LIBOR rates and the barrier is given as an implicit surface in the LIBOR co- 
ordinates. We also note that there is the space domain G in the phase space of 
the SDEs ()6.3|) corresponding to the interval (0, R up ) on the swap-rate semi- line. 
As usual, the corresponding space-time domain Q := [0, To) x G. 

For test purposes, let us introduce an analytical approximation for the barrier 
swaption. To this end, we note that under the Swap Market Model (SMM, see 
details in (4J [18j [21]) the barrier swaption pricing problem admits the closed- 
form solution (cf. (|4.3jl ) 

N 

V swapti on(0) = 6j2 P (°, T i) {Rswap{0) [${6 + {R swap {0) / K, v Rtaap )) (6.4) 
-<P(5+(R swap {0)/R up ,v Rsuiap ))] 

-K [$(5-(R swap (0)/K, v Rsw J) - <5>{5-(R swap (0)/R upi v R _ ap ))] 
-H [HS + (R 2 up /(KR swap (0)),v Rsw J) - H8 + (R up /R swap (0),v Raw J)] 
{0)<P{5-(R up 2 /{KR 

swap 

-H6-(Rup/Rswap{0),V Rsw J)]}, 

where 6± are from (|4.4|) . 

and u Rsmap (s) is the instantaneous volatility of the log-normal dynamics of the 
swap rate. 

Using Rebonato's formula [H p. 283], we can approximately compute the 
"approximate" volatility v R M for the LMM analogous to the volatility v Rs 
in the SMM entering (|B~4|) as 



^ ^(0)^(0)^(0)^(0)^. r T ° 

2^ TZ 77^72 L / o-i{s)(T j (s)ds, (6.5) 



. L MM _ ? 

/■./—ii (Rswap (0)) 
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where 

1-1/ n (i+^(o)) 
5^i/n(i+»(o)) 

k=0 3=° 

The quantity vj^ MM can be used as a proxy for Ur„ in (|6.4I) to compute 
approximated barrier swaption prices under LMM. We will check in our nu- 
merical experiments whether an approximation obtained by our algorithm is 
consistent with this analytical approximation. 

6.1 Algorithm 

Here we exploit Algorithm 12.11 We choose a time step h > so that M = Tq/H 
is an integer. Again inside the domain G we use the weak Euler scheme to 
simulate trajectories of the log LIBOR rates (|6.3I) : 

' SL j 

ln4 +1 = \nL i h + a i (t k )Y,r—*7jPii<rj(tk)h (6-6) 

J=0 1 + 6L k 



I 2 

-o(°"i(*fc)) h + ai(t k )Vh £ U itj (; jk+1 , 

1 3=0 

i = 0,...,N-1, 



where £j k are mutually independent random variables distributed by the law 
P(£ = ±l) = l/2. 



For a fixed t k , we denote by In L k the point with coordinates In L k , In L I 



lnL^ -1 , i.e. lnL k = (lnL k ,lnL\, . . . ,lnL^ _1 ) T . As before, we follow the ran- 
dom walk constructed by (|6.6[) until we reach the boundary zone St kj h- Algorith- 
mically, it implies that we implement a check at each step whether the current 
position of the random walk is in the boundary zone S tk .h- More precisely, we 
evaluate at time t k whether the current position In L k is such that the maximum 
increment from point In L k according to all possible realizations of (|6.6p at the 
next time level t k +i results in the state of the random walk below the barrier, 
i.e. in the domain G. 
Introduce 

InifeMox = max In i*. 

i 

and 

lnife+i = In Lk,Max + ° 2 MaxhN - ^cr 2 M ax h + VMaxVhN, (6.7) 
where (JMax = max^fe cr,(ife). Using the fact that 

Rswap{L k +l, Lk+l) = L k+ i, 
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one can see that the current position of the random walk InLfc is inside the 
domain G if the following condition is satisfied 



Algorithmically, we do the following. If condition (|6.8p is true, we evaluate the 
next position of the random walk at tf.+i according to (|6.6[) and continue further 
with the algorithm unless k + 1 = M (i.e., we have reached the maturity time 
To of the swaption). 

We note the condition (I6.8[) is computationally easy to evaluate but it is 
rather rough. Once this condition fails, we check a finer but computationally 
more expensive condition based on the maximum increments from each of L l (t k ) 
towards the boundary: 



R 



swap 



+ <Jo(t k )(TMax,kh + a (t k )VNh),Ll(l + 2a 1 (t k )a Max ,kh 



+a 1 {t k ) y /(N-l)h),...,L 



N-l 



(1 + N(J j\r-l(tfc)<7 Max,kh + (J N-l {t k )Vh)) < R up , 

(6.9) 



where (JMax,k = max^ aj(t k ). If condition (|6.9[) holds, we again carry on to the 
next time step t k+ i using (I6.6[) and continue further with the algorithm unless 
k + 1 = M. 

If both conditions (I6.8P and (|6.9[) fail, the random walk has reached the 
boundary zone St k ,hi where as before we apply the different procedure which re- 
quire us to find the projection lni^ := (LZ'°, LZ' 1 , ■ ■ ■ , L'^' N ~ 1 ) J of the current 
position In L k on the boundary given as the implicit function of the spanning 
LIBOR rates: 



lnR s 



In i?„ 



(6.10) 



For completeness of the exposition, let us discuss how the projection In L\ 
can be simulated before we return to the description of the algorithm. The 
problem of finding point In L£ is equivalent to finding the minimum value of the 
function 



InLi. I 2 



InL? 



lnL° fc 



r+-+( 



N-l 



InL 



N-l 



(6.11) 



subject to the constraint 



/ 



hi 



N-l , , \ 

n(i+«: j )-i 

J=0 v ' 



N-1 



N-l 



In R y 



(6.12) 



We regard hiLT' 1 , 



InL 71 "' as independent variable in the constraint equa- 
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tion (|6 . 1 2[) and write In L k 
( 



7T.0 



as 



InL? 



In 



R 



up 




\ 



N-l 

n 



(6.13) 



/ 



Hence the minimization problem is reduced to finding the point In L^' 1 , . . . , In L^' N 1 
at which the function | lnL£ — ln.L/.| 2 from (|6.11l) with lnL^'° from (|6.13l) has 
its minimum value. This optimization problem can be solved using standard 
procedures, e.g. the MATLAB function "lsqnonlin()" . 

Let us now continue with the description of the algorithm. When In Lk is in 
the boundary zone, we either stop the chain at In L£ with probability p : 



|lnL£-lnL fe | 



where 



\Vh = VN(a 2 Max hN--a 2 Max 1 



O Ma 



In LV In Lb 



I In L7 — In Lk 



(6.14) 



(6.15) 



with 



or we jump inside the domain Q to the point In Lk + \y/h-r 

probability 1 — p, apply the Euler step (|6.6|) to evaluate In Lk+i and continue 
further with the algorithm unless k + 1 = M. 

The outcome of simulating each trajectory is the point {t >Cl lnL^). In Algo- 
rithm [63] we present the pseudocode for simulating a single trajectory based on 
the algorithm we described above. 



6.2 Numerical results 

We give some results for pricing a barrier swaption by Algorithm 16.11 We 
consider the barrier swaption with the initial LIBOR curve flat at 5%, constant 
volatility ai(t) at 10% and the following parameters: T = 10, T* = 20, K — 
0.01, R up = 0.075, 6= 1, fi = 0.1. The simulations use 10 6 Monte Carlo runs. 

The pricing problem for the barrier swaption (I6.ip does not admit a closed- 
form solution. We used the barrier swaption price found with Algorithm 16.11 
with h = 0.01 and 10 7 of Monte Carlo runs as the reference solution. This 
reference price is 0.15506 with the Monte Carlo error 1.42 x 10" 4 , which gives 
half of the size of the confidence interval for the corresponding estimator with 
probability 0.95. The analytical approximation based on (|6.4I) and (16.5[) yields 
the price of the barrier swaption 0.15556. 

We present results of the experiments in Table loTTl As in the previous sec- 
tion, the error column values before "±" are estimates of the bias computed 
using the reference price value and the values after "±" reflect the Monte Carlo 
error with probability 0.95. The "mean exit time" is the average time for ap- 
proximate trajectories to exit the space-time domain Q. It is clear that the 
results demonstrate the expected first order of convergence. 
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Algorithm 6.1 Barrier swaption: Pseudocode for simulating a single trajectory 
FOR k = 1 to M 

calculate InL^+i by (|6.7|) 

IF (|6.8[) is true: calculate InLfe+i by (|6.6p 

ELSE 

IF (|6.9[) is true: calculate In Lk+i by (|6.6p 
ELSE 

solve the minimisation problem (|6. 11 [) with lnL^'° from (|6.13[) 
generate u ~ Unif[0, 1] 
calculate probability p by (|6.14[) 
IFu<p 
break 

ELSE > 

SET In L k to In L fe + XVh ^^ '"^ 

| In LJ — In in, | 

calculate In Lfe+i by (|6.6p 
ENDIF 
ENDIF 
ENDIF 
ENDFOR 



Table 6.1: Performance of Algorithm 6.1 for the barrier swaption. 



h 




error 




mean exit time 


0.25 


1.01 x 10 


~ 2 ±4.33 x 10 


-4 


9.36 


0.2 


8.08 x 10 


" 3 ±4.37 x 10 


-4 


9.40 


0.125 


5.15 x 10 


- a ±4.42 x 10 


-4 


9.46 


0.1 


4.15 x 10 


" 3 ±4.44 x 10" 


-4 


9.48 


0.0625 


2.58 x 10 


" 3 ±4.47 x 10 


-4 


9.51 


0.03125 


1.03 x 10 


" a ±4.49 x 10 


-4 


9.54 
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